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Abstract — In this paper, numerical methods are suggested to 
compute the discrete and the continuous spectrum of a signal 
with respect to the Zakharov-Shabat system, a Lax operator un- 
derlying numerous integrable communication channels including 
the nonlinear Schrodinger channel, modelling pulse propagation 
in optical fibers. These methods are subsequently tested and their 
ability to estimate the spectrum are compared against each other. 
These methods are used to compute the spectrum of various 
pulses commonly used in the optical fiber communications. It is 
found that the layer-peeling and the spectral methods are suitable 
schemes to estimate the nonlinear spectra with good accuracy. 
To illustrate the structure of the spectrum, the locus of the 
eigenvalues is determined under amplitude and phase modulation 
in a number of examples. It is observed that in some cases, 
as signal parameters vary, eigenvalues collide and change their 
course of motion. The real axis is typically the place from which 
new eigenvalues originate or are absorbed into after travelling a 
trajectory in the complex plane. 

Index Terms — Fiber-optic communications, forward nonlinear 
Fourier transform, Zakharov-Shabat spectral problem, numeri- 
cal methods. 

I. Introduction 

I HE nonlinear Fourier transform (NFT) of a signal q(t) is 
a pair of functions: the continuous spectrum q(X), Ael, 
and the discrete spectrum q(Xj), ^Xj > 0, j — 1,...,N. 
The NFT arises in the study of integrable waveform channels 
as defined in Part I [lj. In such channels, signals propagate 
(in a potentially complicated manner) according to a given 
integrable evolution equation, whereas the nonlinear Fourier 
transform of the signal propagates according to a (simple) 
multiplication operator. 

In [Part I], we proposed nonlinear frequency-division mul- 
tiplexing (NFDM), a scheme that uses the nonlinear Fourier 
transform for data communication over integrable channels. 
NFDM extends traditional orthogonal frequency division mul- 
tiplexing (OFDM) to channels generatable by a Lax pair. An 
example is the optical fiber channel, where signal propagation 
is modelled by the (integrable) nonlinear Schrodinger (NLS) 
equation. In general, the channel input-output relations in the 
NFT domain are (see [Part I]) 

Y(X)=H(X)X(X) + Z(X), 
Y(X j ) = H(X j )X(X j ) + Z(X j ), 
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where X(X) and X(Xj) are continuous and discrete spectra at 
the input of the channel, Y(X) and Y(Xj) are spectra at the 
output of the channel, and Z(X) and Z(Xj) represent noise. 
The channel filter H(X) for the NLS equation is given by 
H(X) = exp(-4?A 2 z). 

NFDM is able to deal directly with nonlinearity and dis- 
persion, without the need for additional compensation at 
the transmitter or receiver. In this scheme, information is 
encoded in the nonlinear spectrum at the channel input, and 
the corresponding time-domain signal is transmitted. At the 
receiver, the NFT of the received signal is computed, and the 
resulting spectra Y(X) and Y(Xj) are subsequently used to 
recover the transmitted information. 

Similar to the ordinary Fourier transform, while the NFT can 
be computed analytically in a few cases, in general, numerical 
methods are required. Such methods must be robust, reliable 
and fast enough to be implemented in real time at the receiver. 
In this paper, we suggest and evaluate the performance of a 
number of numerical algorithms for computing the forward 
NFT of a given signal. Using these algorithms, we then 
perform extensive numerical simulations to understand the 
behavior of the nonlinear spectrum for various pulse shapes 
and parameters commonly used in data communications. 

We are aware of no published work presenting the NFT 
of various signals numerically, for many pulse shapes and 
parameters. Such work is necessary to clarify the structure 
of the nonlinear spectrum and help in its understanding. In 
part, this has been due to the fact that the NFT has largely 
remained a theoretical artifice, and practical implementation 
of the NFT as an applied tool has not yet been pursued in 
engineering. 

We review the relevant literature in Section HI] and suggest 
new schemes for the numerical evaluation of the NFT. Al- 
though these methods are general and work for the AKNS 
system for the purpose of illustration, we specialize 
the AKNS system to the Zakharov-Shabat system. All these 
methods are put to test in cases where analytical formulae exist 
and are compared with one another in Section [V] Only some 
of these methods will be chosen for the subsequent numerical 
simulations; these are the layer-peeling method, Ablowitz- 
Ladik integrable discretization, and the spectral matrix eigen- 
value scheme. These methods are used in the next sections 
to numerically compute the nonlinear Fourier transform of 
a variety of practical pulse shapes encountered in the data 
communications. 
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II. The Nonlinear Fourier Transform 

Details of the nonlinear Fourier transform can be found in 
[Part I]. Here we briefly recall a few essential ingredients re- 
quired in the numerical computation of the forward transform. 
As noted earlier, we illustrate numerical methods in the context 
of the Zakharov-Shabat system, which is a Lax operator for 
the nonlinear Schrodinger equation. 

For later use, we recall that the slowly-varying complex 
envelope q(t, z) of a narrow-band small-amplitude signal 
propagating in a dispersive weakly-nonlinear medium, such 
as an optical fiber, satisfies the cubic nonlinear Schrodinger 
equation. By proper scaling, the equation can be normalized 
to the following dimensionless form in 1 + 1 dimensions: 



qtt+2\q\ 2 q. 



(1) 



Here t denotes retarded time, and z is distance. 

The NFT for an integrable evolution equation starts by 
finding a Lax pair of operators L and M such that the 
evolution equation arises as the compatibility condition L z = 
[M, L] = ML - LM. For the NLS equation, we may take 
operator L as 



L = J 



(2) 



(The corresponding M operator can be found in [Part I].) 

The NFT is defined via the spectral analysis of the L 
operator, given in this paper by d2). The spectrum of L is 
found by solving the eigenproblem Lv = Xv, where A is 
an eigenvalue of L and v is its associated eigenvector. It 
can be shown that the operator L in (O has the isospectral 
flow property, i.e., its spectrum is invariant even as q evolves 
according to the NLS equation. 

The eigenproblem Lv = Xv can be simplified to 

Note that the z-dependence of q is suppressed in (01 (and 
throughout this paper), as this variable comes into play only 
in the propagation of the signal, not in the definition and 
computation of the NFT. 

Assumption 1. Throughout this paper we assume that (a) 
q G L 1 (M), and (b) q(t) is supported in the finite interval 
[T X ,T 2 ). ' □ 

The set of eigenvectors v associated with eigenvalue A in 
@ is a two-dimensional subspace E\ of the continuously 
differentiable functions. We define the adjoint of a vector 
v = [vi(t),v 2 (t)] T as v = [vUt),-vl(t)} T . If v(t,X*) is 
an element of Ex*, then v(t, A*) is an element of E\. It can 
be shown that any pair of eigenvectors v(t,X) and v(t,X*) 
form a basis for E\ [Part I]. Using Assumption 02b), we can 
select an eigenvector A) to be a solution of (01 with the 
boundary condition 



v\T 2 ,X) 



The basis eigenvectors v 1 and v 1 are called canonical eigen- 
vectors. 



Having identified a basis for the subspace E\, we can 
project any other eigenvector v 2 on this basis according to 



v 2 (t, A) = a(A)w 1 (i, A) + biXyit, A). 



(4) 



Following Assumption |TJb), a particular choice for v 2 is made 
by solving the system 



"jA q(t) 



q*(t) jX V > w(Tl ' A) 



(5) 



in which we dropped the superscript 2 in v 2 for convenience. 
By solving (0 in the interval [Ti,T2] for a given A and 
obtaining v(T 2 , A), the nonlinear Fourier coefficients a(A) and 
b(X) can be obtained by considering (@J at t = T 2 . The 
resulting coefficients obtained in this manner are 



a(\)=v 1 (T 2 )e> XT \ 
b{X) = v 2 (T 2 )e 



-jAT 2 



(6) 



The NFT of a signal q(t) consists of a continuous spectral 
function defined on the real axis Ael 



9(A) 



6(A) 
a(XY 



X e 



and a discrete spectral function defined on the upper half 
complex plane C + = {A : 3(A) > 0} 

HA;) 



da(A)/dA| A=A , 



, 3 = !>••• ,N, 



where Xj are eigenvalues and correspond to the (isolated) zeros 
of a(A) in C + , i.e. a(Xj) = 0. 

From the discussions made, in order to compute the nonlin- 
ear spectrum of q(t), the system of differential equations (0 
needs to be solved in the interval [Ti,T 2 ]. Except for special 
cases, © needs to be solved numerically. 

Numerical methods for the calculation of the forward non- 
linear Fourier transform are divided into two classes in this 
article: 

1) Methods which estimate the continuous spectrum by 
directly integrating the Zakharov-Shabat system; see 
Section |III] 

2) Methods which find the (discrete) eigenvalues. Two 
approaches are suggested in this paper for this purpose. 
Similar to the continuous spectrum estimation, we can 
integrate the Zakharov-Shabat system numerically and 
obtain a(A). To find zeros of a(A), the scheme is often 
supplemented with a search method to locate eigenvalues 
in the upper half complex plane. One can also discretize 
and rewrite the Zakharov-Shabat system in the interval 
[Ti , T 2 ] as a (large) matrix eigenvalue problem; see 
Section [IV] 

We begin by discussing methods which estimate the con- 
tinuous spectrum. 

III. Numerical Methods for Computing the 
Continuous Spectrum 

In this section, we assume that A e M is given and provide 
algorithms for calculating the nonlinear Fourier coefficients 
a(A) and b(X). The continuous spectral function is then easily 
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computed as the ratio q(X) = b(\)/a(X), This process can be 
repeated to compute the spectral amplitudes for any desired 
finite set of continuous frequencies A. 

A. Forward and Central Discretizations 

The most obvious method to attempt to solve (0 is the 
first-order Euler method or one of its variations []3]- 

Recall that the signal q(t) is supported in the finite time in- 
terval [TijTa], and partition this interval uniformly according 
to the mesh T x < Ti + e < ■ ■ ■ < T x + Ne = T 2 with size N, 
i.e., with e = (T 2 - T x )/N. Let q[k] = q(ke) and let 

Integrating both sides of (O from ke to (fc + l)e and assuming 
that the right hand side is constant over this interval, we get 

v[k + 1] = v[k] + eP[k]v[k], k = 0,...,N, 

»[0]=(j)e-^. (8) 

Equation ^ is iterated from k = to k = N to find v[N}. 
The resulting vector is subsequently substituted in (|6) to obtain 
a(A) and 6(A). 

We have implemented the Euler method for the calculation 
of the nonlinear Fourier transform of a number of pulse shapes. 
Unfortunately, the one-step Euler method does not produce 
satisfactory results for affordable small step sizes e. 

One can improve upon the basic Euler method by consid- 
ering the central difference iteration J3J, 

v[k + 1] = v[k - 1] + 2eP[k]v[k). 

This makes the discretization second-order, i.e., the error 
v(Ti + ke) — v[k] is of order 0(e 2 ). Here an additional initial 
condition is required too, which can be obtained, e.g., by 
performing one step of the regular forward difference (0. 

B. Fourth-order Runge-Kutta Method 

One can also employ higher-order integration schemes such 
as the Runge-Kutta methods. Improved results are obtained 
using the fourth-order Runge-Kutta method |4j, 0. How- 
ever it takes significant time to estimate the spectrum using 
such higher order numerical methods in real-time. Since the 
method, with its typical parameters, is quite slow and does not 
outperform some of the schemes suggested in the following 
sections, we do not elaborate on this method here; see [3 1 for 
details. However, for comparison purposes we will include this 
scheme in our numerical simulations given in Section [V] 

C. Layer-peeling Method 

In Section IV. C of [Part I], we have calculated the nonlinear 
spectra of a rectangular pulse. One can approximate q(t) as a 
piece-wise constant signal and use the layer-peeling property 
of the nonlinear Fourier transform to estimate the spectrum of 
any given signal. Let a[k] and b[k] be the nonlinear Fourier 
coefficients of q(t) in the interval [T x ,ke), and x[k] and y[k] 



coefficients in the small (rectangular) region [ke, ke + e). The 
iterations of the layer-peeling method read 

(a[k + l],b[k + l]) = {a[k],b[k}) o (x[k],y[k]), (9) 
(o[0],6[Q]) = (l,0), 

where the o operation is defined as in [1| 

a[k + 1] = a[k]x[k] - b[k]y[k], 
b[k + l] = a[k]y[k] + b[k]x[k], 

in which 

x[k] = fcos(Ue) - j-^sin(De)^ e JA(t[fe]-t[fc-i]) ) 

y[k] = ZA s ^ De ) e -3Mt[k]+t[k-l])^ 

an d x[k}(\) = x*[k}(\*), y[k](X) = y*[k](X*), D = 
y/X 2 + \q[k]\ 2 . The desired coefficients are obtained as a := 
a[N] and b := b[N]. Note that the exponential factors in x[k] 
and y[k] enter ( TTOb in a telescopic manner. As a result, for the 
numerical implementation, it is faster to drop these factors 
and just scale the resulting a[N] and b[N] coefficients by 
exp(jA(T 2 -Ti)) and exp(-jA(T 2 +Ti)), respectively. This, 
however, reduces the accuracy as it involves the product of 
large and small numbers. 

We are motivated by (6) in which the layer-peeling identity 
(O is mentioned as a property of the nonlinear Fourier 
transform. An equivalent presentation of this method is given 
in Q as well. 

Note further that a different numerical method, but with 
the same name (layer-peeling), exists in geophysics and fiber 
Bragg design [8|; however this method is not related to the 
problem considered here. 

We shall see in Section [V] that the layer-peeling method 
gives remarkably accurate results in estimating the nonlinear 
Fourier transform. 



D. Crank-Nicolson Method 

In the Crank-Nicolson method, the derivative of the evolu- 
tion parameter is approximated by a finite-difference approx- 
imation, e.g., forward discretization, and other functions are 
discretized by taking their average over the end points of the 
discretization interval: 

V[k + 1] e ~ V[k] = \ (P[k]v[k] + P[k + l]v{k + 1]) , 

where P[k] is defined in (0. This implicit iteration can be 
made explicit 

v[k + 1] = (I - e -P[k + l])" 1 ^ + ^P[k])v[k], 
k = 0,...,N 

with initial condition ((HJ. As we will see, this simple scheme 
too gives good results in estimating the nonlinear spectrum. 
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E. The Ablowitz-Ladik Discretization 

Ablowitz-Ladik discretization is an integrable discretization 
of the NLS equation in time domain [9|. In this section, we 
suggest using the Lax pairs of the Ablowitz-Ladik discretiza- 
tion of the NLS equation for solving the Zakharov-Shabat 
eigenproblem in the spectral domain. 

Discretization sometimes breaks symmetries, making the 
discrete version of an integrable equation no longer integrable. 
A consequence of symmetry-breaking is that quantities which 
are conserved in the continuous model may no longer be 
invariant in the discretized equation. A completely integrable 
Hamiltonian system with an infinite number of conserved 
quantities might have a discretized version with no, or few, 
conserved quantities. The discrete equation therefore does not 
quite mimic the essential features of the original equation if 
the step size is not small enough. 

However, for some integrable equations, discretizations ex- 
ist which are themselves completely integrable Hamiltonian 
systems, i.e., they possess an infinite number of conserved 
quantities and are linearizable by a Lax pair, and therefore 
are solvable by the nonlinear Fourier transform. Such de- 
velopments exist for the NLS and Korteweg-de Vries (KdV) 
equations. 

For the NLS equation, the integrable discrete version was 
introduced by Ablowitz and Ladik [9 |. To illustrate the general 
idea, let us replace l±j Xe for small e with e ±jAe in the forward 
discretization method <[8J (the opposite of what is usually done 
in practice). Let z — e~ jAe represent the discrete eigenvalue, 
Q[k] = q[k]e, and 



R[k] = 



-Q* 



Q[k]\ 



The Ablowitz-Ladik iteration is 



v[k + 1] = R[k]v[k), v[0} = 



(10) 



(ID 



Under the z transformation, the upper half complex plane in 
A domain is mapped to the exterior of the unit circle in the 
z domain. The continuous spectrum therefore lies on the unit 
circle \z\ = 1, while the discrete spectrum lies outside of the 
unit circle \z\ > 1, 

One can rewrite the i?-equation ( TTOb in the eigenvalue form 
Lv[k] = zv[k], with the following L operator 



L 



Z 

-Q*[k-1] 



-Q[k] 
a[k - \\Z~ X 



(12) 



where a[k] = 1 + |Q[fc]| 2 , and Z is the shift operator, i.e., 
Z~ 1 x[k] — x[k — 1], k € Z. To the first order in e, a[k] w 1 
and (fT2l) can be simplified to 



L 



Z 
-Q*[k] 



-Q[k}\ 



1 ■ 



(13) 



Given the L operator ([T3V one can consider the M operator 
of the continuous NLS equation and modify its elements such 
that the compatibility equation L z = [M, L] represents a 
discretized version of the NLS equation. It is not hard to verify 



that after doing so we are led to an M operator resulting in 
the following discrete integrable NLS equation 



.dq[k] _ q[k + 1] - 2q[k] + q[k - 1] 



3- 



dz 



+ \q[k]\ 2 (q[k + l}+ q [k-l}). 



(14) 



Here the space derivative remains intact and the signal q[k] 
is discretized in time, in such a way that the nonlinearity 
is somehow averaged among three time samples. In the 
continuum limit e — > 0, ( fT4l i approaches the continuous NLS 
equation and its merits lie in the fact that it is integrable for any 
e, not just in the limit e — > 0. For example, soliton pulses can 
be observed in this model for any e. The equation has its own 
infinite number constants of motion, approaching integrals of 
motion in the continuum limit. The operator M which leads 
to ( fT4t . and the details of the nonlinear Fourier transform for 
CE!) can be found in l9l. ifTUl. 

We conclude that the Ablowitz-Ladik discretization can be 
used not only as a means to discretize the NLS equation in 
the time domain, but also as a means to solve the continuous- 
time Zakharov-Shabat system in the spectral domain. This is 
a non-finite-difference discretization, capable of dealing with 
oscillations exp(±j Ai) in the Zakharov-Shabat system, which 
greatly enhances the accuracy of the one-step finite-difference 
methods. 

Following the Tao and Thiele's approach J6] and J4), we 
can also normalize the R[K] matrix 



v[k + l] 



1 



v/i + |QWI 2 V-OM 



Q[k}\ 
z- 1 



v[k]. (15) 



The scale factor does not change the spectrum significantly, 
since it is cancelled out in the ratios q = b/a and q = b/a', and 
also its effects are second order in e. However, numerically, 
normalization may help in reducing the numerical error. In 
subsequent sections, we refer to ([Til as the Ablowitz-Ladik 
method (AL1) and to (T5[ as the modified Ablowitz-Ladik 
method (AL2). 

IV. Methods for Calculating the Discrete 
Spectrum 

In order to compute the discrete spectrum, the zeros of 
a(A) in the upper half complex plane must be found. One 
way to visualize this is to assume a two-dimensional mesh in 
C + and determine a at all mesh points. Discrete eigenvalues 
are then easily identified by looking at the graph of |a(A)|; 
in many cases they correspond to deep and narrow "wells" 
corresponding to the zeros of the magnitude of a. 

As noted earlier, two types of methods are suggested to 
calculate the point spectrum. 

1) One can use the integration-based algorithms mentioned 
in Section [III] which calculate nonlinear Fourier coeffi- 
cients, and search for eigenvalues using a root finding 
method, such as the Newton-Raphson method. Such 
methods require good initial points and one needs to 
be careful about convergence 

2) It is also possible to rewrite the spectral problem for 
an operator as a (large) matrix eigenvalue problem. The 
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point spectrum of the operator can be found in this way 
too. 



A. Search Methods 

To calculate the discrete spectral amplitude q = b/a', we 
require da/dX as well. As we will show, information about 
the derivative of a can be updated recursively along with 
the information about a, without resorting to approximate 
numerical differentiation. 

Recall that the nonlinear Fourier coefficient a(A) is given 
by © 

o(A) =v 1 [N}e lXt[N] . 
Taking the derivative with respect to A, we obtain 

^ = {{vim > +j T 2Vl [N})e^W. 

We can update the derivative information dv/dX along with v. 
In methods of Section [TTT1 the transformation of eigenvectors 
from t[k] to t[k + 1] can be generally represented as 

v[k + 1] = A[fc]u[fc], 

for some suitable one-step update matrix Ak (which varies 
from method to method). Differentiating with respect to A and 
augmenting v with v' = dv/dX, we get the iterations 



v[k + 1] = A[fc]v[fe], 
v'[k + 1] = A'[k)v[k) + A[k]v'[k], 

with initial conditions 



(16) 
(17) 



-j\t[0] 



-M0]\ e _ jXt 



[0] 



The derivative matrix A' depends on the method used. 
For the forward discretization scheme: 



A' = Mi = 



-3 
j 



For the Crank-Nicolson scheme: 



A' = ^(I + M 2 ) 1 {l+{I + M 2 ) 1 {I-M 2 ))M l , 



where 

\]X -\q[k] 
}q*[k] -\jX 

For the Ablowitz-Ladik method: 



M, 



A 1 = 



~3Z 







jz- 

The desired coefficients are obtained at k = N as follows 

a(A) = Vl [N]e jXt W, 
6(A) =v 2 [N]e- jXt W, 
a'(A) - (v{[N] + j*[JV]«i[JV]) e* x *W . 



Similarly, the layer-peeling iteration can be augmented to 
update a' (A) as well: 



a[k+l] -- 
b[k + 1] = 
a'[k + l] -- 
b'[k+l] - 
a[0]-- 

where 

x[k] = 
y[k] = 

x'[k] = 
y'[k] = 



a[k}x[k] - b[k]y[k], 
a[k]y[k] + b[k]x[k], 

a'[k}x[k] + a[k]x'[k] - (b'[k]y[k] + b[k]y'[k]), 
a'[k}y[k] + a[k}y'[k] + b'[k]x[k] + b[k]x'[k], 
1, 6[0] = o'[0] = 6'[Q] = 0, (18) 



cos(De) - j— sm(De), 

--^-cospe) - — - ) sm(De), 
-J^ (Decos(De)-sm(De)) 



The expressions for x' [k] and y' [k] are similar, with j replaced 
with — j and q[k] replaced with — q*[k]. 

With the derivative information being available, the Newton- 
Raphson method is a good scheme to search for the location 
of the (discrete) eigenvalues. The iteration for the complex- 
valued Newton-Raphson scheme is 



Afe+i — Xk — «fc 



a(Afc) 
a'(Afe) 



(19) 



where ak is some step size modifier; usually ak = 1. The 
iteration stops if A& is almost stationary, i.e., if \a^r\ < 5 for 
a small 5. In practice, the quadratic convergence of the scheme 
is often very fast and occurs in just a few iterations. 

In data communications, since noise is usually small, the 
points in the transmitted constellation can (repeatedly) serve 
as the initial conditions for (19[ . In this case, convergence 
is usually achieved in a couple of iterations. For an unknown 
signal, random initial conditions are chosen. In either case, one 
or more sequence of Newton iterations have to be performed 
for any single eigenvalue. 

To make sure that all of the eigenvalues are found, we can 
check the trace formula for n — 1,2,3 [Part I]. The trace 
formula is a time frequency identity relating the hierarchy of 
infinitely many conserved quantities to the spectral compo- 
nents. 

In general, the trace formula represents a time domain 
conserved quantity as the sum of discrete and continuous 
spectral terms: 



E {k) 



rpi k ) , EA 
-^disc ' -^conti 



where 



E, 



(k) 



A fc - 1 log(l + | ( ?(A)| 2 ) dX, 
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and where are time domain conserved quantities (func- 
tionals of the signal). The first few conserved quantities are 
energy 



\q(t,z)\ 2 dt, 



momentum 



oo 

E = 2] J q{t) — dt > 



and Hamiltonian 



B (3) 




dq(t) 



dt 



dt. 



For n = 1, the trace formula is a kind of Parseval's identity 
(Plancherel's theorem), representing the total energy of the 
signal in time as the sum of the energy of the discrete and 
continuous spectral functions. When satisfied, the Parseval's 
identity ensures that all of the signal energy has been ac- 
counted for. 

We first calculate the continuous spectrum and its "energy 
terms" for a sufficiently fine mesh on the real A axis. The 
energy difference E en:ol — \\E^ — -Etontll gives an estimate 
on the number of missing eigenvalues. When a new eigenvalue 
A is found, this error is updated as E errol := \\E^ — E^t ~~ 
r^A ||. The process is repeated until E aTOI is less than a small 
prescribed tolerance value. 

To summarize, given the signal q(t), its nonlinear Fourier 
transform can be computed based on Algorithm 1 . 



B. Discrete Spectrum as a Matrix Eigenvalue Problem 

The methods mentioned in Section IIV-AI find the discrete 
spectrum by searching for eigenvalues in the upper half com- 
plex plane. Sometimes it is desirable to have all eigenvalues 
at once, which can be done by solving a matrix eigenvalue 
problem [5]. These schemes obviously estimate only (discrete) 
eigenvalues and do not give information on the rest of the 
spectrum. Since the matrix eigenvalue problem can be solved 
quickly for small-sized problems, it might take less computa- 
tional effort to compute the discrete spectrum in this way. In 
addition, one does not already need the continuous spectrum to 
estimate the size of the discrete spectrum. On the other hand, 
for large matrices that arise when a large number of signal 
samples are used, the matrix eigenproblem (which is usually 
not Hermitian) is slow and it may be better to find the discrete 
spectrum using the search-based methods. The matrix-based 
methods also have the disadvantages that they can generate a 
large number of spurious eigenvalues, and one may not be able 
to restrict the algorithm for finding eigenvalues of a matrix to 
a certain region of the complex plane. 

Below we rewrite some of the methods mentioned in the 
Section [HI] as a regular matrix eigenvalue problem, for the 
computation of the discrete spectrum. 



Algorithm 1 Numerical Nonlinear Spectrum Estimation 
Sample the signal at a sufficiently small step size e. 
Fix a sufficiently fine mesh M on the real A axis, 
for each A e M do 

Iterate (O from k — to k = N to obtain v[N]. 

Compute the continuous spectral amplitude 

vi[N\ 

end for 

Initialize the error e = \\E — _E cont 1 1 . 
while |e| > ei do 

Choose Ao 6 T> randomly, where I? is a prescribed region 

inC+. 

Set i = 0; 

repeat 

Iterate (|T6j-([T7j from k = to k = N to obtain v[N] 
and u'[iV], and perform a Newton-Raphson update 

vi[N] 



A 



Xi - AA, AA 



a. 



v[[N} + jt[N} Vl [N} 



if If \ i+ i£V then 

choose Ao € T> randomly. 

Set i := -1. 
end if 

Set i:=i+l. 
until |AA| < e 2 and i > 

Aj is an eigenvalue and the associated spectral amplitude 
is 

v 2 [N] 



v[[N] + jt[N} Vl [N] 



Update e 



\E — E rn „, — Ea 



-2jX t t[N\ 



where E disc 



[43^,2^,1^]. 
end while 



1 ) Central-difference Eigenproblem: The matrix eigenvalue 
problem can be formulated in the time domain or the frequency 
domain @. Consider the Zakharov-Shabat system in the form 

Lv = Xv 



3 



d_ 

dt 

-q* 



v = Xv. 



(20) 



In the time domain, one can replace the time derivative d/dt 
by the central finite-difference matrix 





/ 





1 





... _1\ 


1 




-1 





1 


... o 


2e 










-1 


1 






1 








-1 o) 



and expand ( f20T > as 
D 

-diag(g*[ 



J 



-diag(#])\ v[k] 



Xv[k\. (21) 



The point spectrum is contained in the eigenvalues of the 
matrix in the left hand side of (fJTJ. 
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Eigenvalues of a real symmetric or Hermitian matrix can be 
found relatively efficiently, owing to the existence of a com- 
plete orthonormal basis and the stability of the eigenvalues. 
In this case, a sequence of unitary similarity transformations 
Ak+i = PAkP T can be designed, using, for instance, the QR 
factorization, the Householder transformation, etc., to obtain 
the eigenvalues rather efficiently IfTTI . 

Unfortunately, most of the useful statements about compu- 
tations using Hermitian matrices cannot usually be generalized 
to non-Hermitian matrices. As a result, the eigenvalues of 
a non-Hermitian matrix (corresponding to a non-self-adjoint 
operator) are markedly difficult to calculate ifPTl . lfl"2l . Run- 
ning a general-purpose eigenvalue calculation routine on (f21~b 
is probably not the most efficient way to get eigenvalues. 
Next we make suggestions to simplify the non-Hermitian 
eigenproblem (l2TT i by exploiting its structure. 

The diagonal matrix in the lower left corner of (12H can 
be made zero by applying elementary row operations and 
using the entries of the D matrix. Since elementary row 
operations, as in Gauss-Jordan elimination, generally change 
the eigenvalues, the corresponding column operations are also 
applied to induce a similarity transformation. In this way, 
an upper-Hessenberg matrix is obtained in O(N) operations, 
as compared to a sequence of Householder transformations 
with 0(N 3 ) operations and 0(N 2 ) memory registers. The 
eigenvalues of the resulting complex upper-Hessenberg matrix 
can subsequently be found using QR iterations. 

2) Ablowitz-Ladik Eigenproblem: We can also rewrite the 
Ablowitz-Ladik discretization as a matrix eigenvalue problem. 
Using the L operator ([T3~l >. we obtain 



-Q*[k 



v\[k + 1] — Q[fc]w 2 [fc] = 
l>i[fc] + a[k - l]v 2 [k - 1] = zv 2 [k], 



(22) 
(23) 



which consequently takes the form 

1]) 



-diag(Q*[fc 



d H Q[fc]) ) v = zv, 



U. 



in which 



m = 



U 2 = 





i 


••• 










1 














1 


V 








V 


( 





a[l] 













a[2] 














\a[N] 






a[N ~ 1] 




/ 



and diag(Q[fc]) = diag (Q[0], • • ■ ,Q[N]), diag(Q*[* - 1]) - 
diag (Q*[N],Q*[0] ■■■ ,Q*[N - 1]). Note that all shifting op- 
erations are cyclic, so that all vector indices k remain in the 
interval < k < N. 

It is not very different to consider the more simplified L 
operator (fT2l instead of ([T3V This corresponds to the above 
discretization with U2 '■= U\ and — Q*[k — 1] replaced with 



— Q*[k). Similarly, one can rewrite the normalized Ablowitz- 
Ladik iteration as a matrix eigenvalue problem. This corre- 
sponds to ( |23l where instead of having a[k — 1], Vi[k + 1] is 
multiplied by a[k] in the first equation, i.e., U\ and U2 are 
interchanged. 

3) Spectral Method: In the frequency domain, one can 
approximate derivatives with the help of the Fourier transform. 
Let us assume that 



a k \ j2k*t 

*) e * , ff(t) 



E 3'2fc7Tt 
Ike t . 



«(*)= E 

where T = T 2 — T\. Then the Zakharov-Shabat system is 



-otk- 



2irk 



2_j Ik-mPm = Aa fe , 



M 



m=-4f 



Thus we obtain 



where a = [a_ m , 



,Q!m] t , (3 = [I3-m,---,(3m] t , n 



■fdiag(-f and 



/ 7o 7-i ■•• 7_m 

71 70 7-1 ••• 7_m 



7*f 7m_i 
2 2 1 

7m 7m_ 



^ 

••• 



7-f +1 7-f 



• • • 7m 7m _i ' ' ' 

2 2 1 

\ • • • 7jv/ 7m _ 1 



7o 



7-i 
7o 



/ 



The point spectrum is thus found by looking at the eigenvalues 
of the matrix 



A 



C. Running Time, Convergence and Stability of the Numerical 
Methods 

The numerical methods discussed in this paper are first- 
order matrix iterations and therefore the running time of all of 
them is O(N) multiplications and additions per eigenvalue. 
This corresponds to a complexity of 0(N 2 ) operations for 
the calculation of the continuous spectrum on a mesh with N 
eigenvalues. The exact number of operations depends on the 
details of the implementation and the memory requirement of 
the method. All iterative methods thus take about the same 
time asymptotically, albeit with different coefficients. 

An important observation is that, while the Fast Fourier 
Transform (FFT) takes 0(A^log 2 A^) operations to calculate 
the spectral amplitudes of a vector with length N at N 
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equispaced frequencies, the complexity of the methods de- 
scribed in this paper to compute the continuous spectrum are 
0(N 2 ). Similarly, it takes 0(NN) operations to calculate the 
discrete spectrum. In other words, so far we do not exploit the 
potentially repetitive operations in our computations. 

It is evident from © that as T 2 — > 00, v\ [k] should grow 
as ~ cxp(— j\T 2 ) so that a(A) is a finite complex number. 
The canonical eigenvector v[k;Ti,T 2 ] thus has an unbounded 
component as T 2 — > 00 (i.e., ||w[fc]|| — > 00). One can, however, 
normalize v\ and v 2 according to 



vie 



u 2 = v 2 e 



and transform (0 to 



Ui 





-q*{t)e 



-2jXt 



q(t)e 




2j\t 



-jXt 



u, u(Ti,A)=(q).(24) 



The desired coefficients are simply a(A) = u\(T 2 ) and 
6(A) = u 2 (T 2 ). Consequently, if one is interested in obtaining 
eigenvectors v[k] in addition to the coefficients a(A) and 6(A), 
the discretization of the normalized system (l24l has better 
stability properties: 

1 Q[k]z 
~Q*[k]z 2k 1 



(a[fc + !],&[&+!]) 



-2k 



(o[0],6[0]) = (l,0). 



(a[k],b[k}) 
(25) 



The nonlinear Fourier coefficients are obtained as a := a[N] 
and 6 := b[N]. The discrete nonlinear Fourier transform 
mentioned in H is thus the forward discretization of the 
normalized Zakharov-Shabat system (|24V 

We are interested in the convergence of v[k] (or a(A) and 
6(A)) as a function of N for fixed values of T\ and T 2 . That 
is to say, we require that the error e = \\v(ke) — v[k]\\ — s- 
as N — > 00 (for fixed T\ and T 2 ). The (global) error in all 
methods described in this paper is at least 0(e), and therefore 
all these methods are convergent. 

Some of these methods are, however, not stable. This is 
essentially because the Zakharov-Shabat system in its orig- 
inal form has unbounded solutions, i.e., ||v(<)|| — >• 00 as 
t —¥ 00. Errors can potentially be amplified by the unstable 
system poles. One should be cautious about the normalized 
system ( f24b as well. For example, forward discretization 
of the normalized system (l25l > gives a first-order iteration 
x[k + 1] = A[fc]x[/c]. The eigenvalues of the matrix A[k] in 
this method are 

Sl , 2 = l±je|#]|. 

It follows that the forward discretization of d24l gives rise to 
eigenvalues outside of the unit disk, |s| > 1. As a result, first- 
order discretization of (l24l are also unstable. In cases where 
\s\ > 1, we can consider normalizing the iterations by dividing 
A[k] by -y/det A[fc] (in the case of d25l l. dividing the right-hand 
side by -Jl + \Q[k]\ 2 ). The resulting iteration has eigenvalues 
inside the unit disk. For e <C 1 the effect is only second-order 
in e, however it helps in managing the numerical error if larger 
values of e are chosen. 

An issue pertinent to numerical methods is chaos. Chaos 
and numerical instability of finite-difference discretizations has 



been observed in [13] for the sine-Gordon equation, which 
is also integrable and shares a number of basic properties 
with the NLS equation. In [14], the authors conclude that 
the standard discretizations of the cubic nonlinear Schrodinger 
equation may lead to spurious numerical behavior. This insta- 
bility is deeply related to the homoclinic orbits of the NLS 
equation, i.e., it occurs if the initial signal q(t, z — 0) is 
chosen to be close to the homoclinic orbit of the equation. 
It disappears only if the step size is made sufficiently small, 
which can be smaller than what is desired in practice. 

It is shown in [14] that the Ablowitz-Ladik discretization 
of the NLS equation (in time) has the desirable property that 
chaos and numerical instability, which are sometimes present 
in finite-difference discretizations of the NLS equation, do 
not appear at all. Though these results are for the original 
time-domain equations, the issue can occur in the spectral 
eigenvalue problem (0 as well, if the signal q[k] is close to 
a certain family of functions (related to s'muit and coswi). 
Therefore, among discretizations studied, the Ablowitz-Ladik 
discretization of the Zakharov-Shabat system is immune to 
chaos and its resulting numerical instability. This is particu- 
larly important in the presence of amplifier noise, where chaos 
can be more problematic. 

V. Testing and Comparing the Numerical Methods 

In this section, we test and compare the ability of the 
suggested numerical schemes to estimate the nonlinear Fourier 
transform (with respect to the Zakharov-Shabat system) of var- 
ious signals. Numerical results are compared against analytical 
formulae, in a few cases where such expressions exist. Our aim 
is to compare the speed and the precision of these schemes 
for various pulse shapes in order to determine which ones are 
best suited for subsequent simulation studies. 

To derive the analytical formulae, recall that the continuous 
spectral function can be written as q(X) — lim y(t, A), in 

t— >oo 

which y(t, A) satisfies HI 

+ q (t)e^Y(t, A) + q*(t)e- 2 ^ = 0, 

y(-oo,A)=0. (26) 
Similarly, one can solve the second-order differential equation 



dt 2 



2jX + - 

q 



dt 



z(-oc,A) = l, ^Z^ =0) (27) 
dt 



and obtain a(A) = lim z(t, A). The zeros of a(A) form the 

t— >oo 

discrete spectrum. 

In the following, the discrete spectrum is found and com- 
pared using the following matrix-based schemes: 

1) central difference method; 

2) spectral method; 

3) Ablowitz-Ladik (AL) discretization (AL1); 

4) AL discretization with normalization (AL2). 

In the matrix-based schemes, the entire point spectrum is 
found at once by solving a matrix eigenvalue problem. 
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The complete spectrum is found using search-based meth- 
ods: 

1) forward discretization method; 

2) fourth-order Runge-Kutta scheme; 

3) layer-peeling methods; 

4) Crank-Nicolson method; 

5) AL discretization; 

6) AL discretization with normalization. 

In search-based methods, the Newton method is used together 
with the trace formula to find both discrete and continuous 
spectra. 

Each of the following pulses is sampled uniformly using a 
total of n samples in a time window containing at least 99.99% 
of the pulse energy. 

A. Satsuma-Yajima pulses 

One signal with known spectrum is the Satsuma-Yajima 
pulse H3 

q(t) = Asech(t). 

Solving the initial value problem (|5]l (or ( |26l i) analytically, the 
following continuous spectral function is obtained |fT31 

r(-j'A + \ + A)T(-j\ + i - A) sin(TrA) 



?(A) 



r 2 (-jA + i) 



cosh(7rA) 



The discrete spectrum is the set of zeros of a(A) i.e., poles 
of q(X) (when analytically extended in C + ). Recalling that 
T(x) has no zeros and is unbounded for x = 0,-1,-2,..., 
it follows that the discrete spectrum consists of N = [A + 
\ — ej eigenvalues, located at (A — (A — . . .. In the 
special case in which A is an integer, A = N, the Satsuma- 
Yajima pulse is a pure N-soliton with N eigenvalues, and the 
continuous spectral function is zero. 

Figs. [T] |2] and [3] give the numerical results for A — 2.7, 
N = 2 10 . 

Fig. Q] shows that, in this example, the spectral and central 
difference methods produce good results among the matrix- 
based methods in estimating the discrete eigenvalues A = 
0.2j, 1.2 j, 2.2j. All methods generate a large number of spu- 
rious eigenvalues along the real axis. This behavior might 
be viewed as a tendency of the algorithms to generate the 
continuous spectrum too. However the spurious eigenvalues do 
not disappear completely even when the continuous spectrum 
is absent (when A is an integer); only their range becomes 
more limited. The spurious eigenvalues across the real axis 
can easily be filtered, since their imaginary part has negligible 
amplitude. The AL methods, with and without normalization, 
produce the same eigenvalues plus another vertical line of 
spurious eigenvalues having a large negative real part. Nor- 
malization in the AL scheme does not make a significant 
difference in this example. 

Fig. |2] shows the accuracy of the various matrix-based 
methods in estimating the smallest and largest eigenvalues 
of q — 2.7sech(i) in terms of the number of the sample 
points N. As the number of sample points N is decreased, the 
spectral and central difference methods maintain reasonable 
precisions, while the accuracy of the AL schemes quickly 
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Fig. 1. Discrete spectrum of the Satsuma-Yajima pulse with A = 2.7 using 
(a) central difference method, (b) spectral method, (c) Ablowitz-Ladik scheme, 
(d) modified Ablowitz-Ladik scheme. 
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Fig. 2. Error in estimating (a) the smallest eigenvalue and (b) the largest 
eigenvalue of the Satsuma-Yajima pulse q(t) = 2.7sech(t) as a function of the 
number of sample points using matrix eigenvalue methods. The Ablowitz- 
Ladik method 1 is the method of Section IIV-B2I with no normalization, and 
the Ablowitz-Ladik method 2 is the same scheme with normalization. 



deteriorates. One can check that in these cases, the error in 
the approximation e jAAt ~ 1 + jXAt becomes large (since 
At > 1). 

The spectral method is generally more accurate than the 
other matrix-based methods. The AL discretizations seem to 
perform well as long as AAi -C 1, i.e., when estimating 
eigenvalues with small size or when N > 200. The AL 
discretization eventually breaks down at about N = 50 as the 
analogy between the continuous and discrete NLS equation 
is no longer justified at such low resolutions, whereas other 
schemes continue to track the eigenvalues to some accuracy. 
In other words, what the AL methods find at such small values 
of N is the spectrum of the discrete soliton-bearing NLS equa- 
tion, which is not a feature of finite-difference discretizations. 
(In fact, it is essential for this algorithm to deviate from the 
finite-difference discretizations as N is reduced, to produce 
appropriate solitons with few samples.) The running time of 
all matrix eigenvalue methods is about the same. 

Search-based methods can be used to estimate the point 
spectrum as well. Here we use the Newton method with 
random initial points to locate eigenvalues in C + . Naturally, 
we limit ourselves to a rectangular region in the complex 
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Fig. 3. Error in estimating the largest eigenvalue of Satsuma-Yajima pulse 
q(t) = 2.7sech(t) as a function of the number of sample points N using 
search-based methods. 



plan, slightly above the real axis to avoid potential spurious 
eigenvalues. Since the number of eigenvalues is not known a 
priori, the continuous spectrum is found first so as to give 
an estimate of the energy of the discrete spectrum. It is 
essential that the continuous spectrum is estimated accurately 
so that a good estimate of the energy of the discrete spectrum 
can be obtained. Once this energy is known, and a suitable 
(rectangular) search region in C + is chosen, the Newton 
method is often able to locate all of the discrete eigenvalues 
using just a few iterations. 

Fig. |3] shows the accuracy of the searched-based meth- 
ods in estimating the largest eigenvalue of the signal 
q(t) = 2.7sech(i). The Runge-Kutta, layer-peeling and Crank- 
Nicolson methods have about the same accuracy, followed 
closely by forward discretization. Since this is the largest 
eigenvalue, the AL schemes are not quite as accurate. As noted 
above, comparison at smaller values of N is not illustrative, as 
the AL estimate quickly deviates from A ma x of the continuous 
signal. 

The Runge-Kutta method, at the accuracies shown in the 
above graphs, is of course very slow, and is not a practical 
method to implement. The running time for the other schemes 
is approximately the same. Search-based methods take an or- 
der of magnitude more time than matrix-based methods when 
N is small. These methods fail when N becomes too small 
(N < 200), since the large error in estimating the energy terms 
of the continuous spectrum negatively influences the stopping 
criteria and consequently degrades the Newton increments. 
For large N, on the other hand, the QR factorization, which 
takes 0(N 3 ) operations in calculating the eigenvalues of a 
matrix, becomes quite slow and restricts the use of matrix- 
based methods. 

The same conclusions are observed for various choices of 
real or complex parameter A. As \A\ is increased, as before, 
the spectral and finite-difference schemes produce the correct 
eigenvalues, and the AL methods generate the same eigenval- 
ues plus an additional vertical strip of spurious eigenvalues. 
The range of the spurious eigenvalues across the real axis 
remains about the same. As the phase of A is increased, 
the true (non-spurious) eigenvalues remain the same in all 
methods (as expected analytically), while some of the vertical 
spurious eigenvalues in the AL schemes move from left to 
right or vice versa. The spectral and finite-difference schemes 
are relatively immune to these additional spurious eigenvalues. 
Normalization of the AL method sometimes produces slightly 



fewer spurious eigenvalues across the real axis, as can be seen 
in Figs. |TJc)-(d). 



B. Rectangular pulse 

Consider the rectangular pulse 

'A iG[T x ,T 2 ] 
otherwise 



(28) 



It can be shown that the continuous spectrum is given by Q] 

D 



A* 

<?(A) - —e 
jA 



2jXt 



1 



.A 



cot(£>(T 2 - T x )) 



where D = ^/\ 2 + \ A\ 2 . To calculate the discrete spectrum, 
the equation ( |27| i is reduced to a simple constant coefficient 
second-order ODE 

^ - 2j\^ + \A\ 2 z = 0, z(Ti) = M'(Ti) = 0. 
It is easy to verify that the eigenvalues are the solutions of 

e 2i(T 2 -Ti)VA3+|Ap = A + V A2 + ijF f2 9) 

\~ \/\ 2 + \A\ 2 ' 

Following the causality and the layer-peeling property of 
the NFT, one can generalize the above result to piece-wise 
constant pulses. This is the basis of the layer-peeling method 
of Section [ITLCl 

Fig. Sta)-(b) show the results of numerically computing 
the discrete spectrum of a rectangular pulse with parameters 
A = 2, T2 = —T\ = 1. The exact eigenvalue is found to 
be A = 1.5713j, by numerically finding the roots of ( |29i > 
using the Newton-Raphson method. No other eigenvalue is 
found under a large number of random initial conditions. 
All methods generate the desired eigenvalue together with 
a large number of spurious eigenvalues across the real axis. 
The central-difference scheme visibly generates fewer spurious 
eigenvalues. The Ablowitz-Ladik schemes produce two more 
eigenvalues with a large negative real part. 

Fig. [5] compares the precision of various methods in esti- 
mating the nonlinear spectrum of the rectangular pulse with 
A = 2, T2 = — Tx = 1. The modified AL scheme performed 
the same as the basic AL scheme, and hence we do not include 
the modified AL scheme in the graphs. 

Convergence to the discrete spectral amplitudes generally 
occurs much more slowly than convergence to the eigenvalues. 
Fig-ISa) shows the precision of various methods in estimating 
the discrete spectral amplitude of the rectangular wave with 
A = 2, T2 = —T\ = 1. It can be seen that convergence does 
not occur until N > 1000. Fig. |6jb) shows the continuous 
spectrum for the same function. All methods produced essen- 
tially the same continuous spectrum, except for some very 
slight variations near zero frequency. 

As \A\ is increased, more eigenvalues appear on the imagi- 
nary axis. The distance between these eigenvalues becomes 
smaller as |A| is increased. All methods produce similar 
results, with the Ablowitz-Ladik methods reproducing the 
purely imaginary eigenvalues at spurious locations with large 
real part. Phase addition has no influence on any of these 
methods, as expected analytically. 
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Fig. 4. Discrete spectrum of the rectangular pulse J28 1 with A = 2, 
T2 = — Ti = 1 using (a) Fourier method, (b) central difference method, 
(c) Ablowitz-Ladik scheme, (d) modified Ablowitz-Ladik scheme. 
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Fig. 5. Error in estimating the largest eigenvalue of the rectangular pulse 
q(t) = 2rect(t) as a function of the number of sample points JV using (a) 
matrix-based methods and (b) search-based methods. 
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Fig. 7. (a) Amplitude profile of a 4-soliton pulse with spectrum <30t . (b) 
Error in estimating the eigenvalue A = 1 + 0.5j. 
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Fig. 8. Discrete spectrum of N-soliton pulse with spectrum 1301 : (a) 
Fourier method, (b) central difference method, (c) Ablowitz-Ladik scheme, 
(d) modified Ablowitz-Ladik scheme. 



C. N-soliton pulses 

We consider an 4-soliton pulse with discrete spectrum 

q(-l + 0.25j) = 1, q(l + 0.25j) = —j, 
g(-l+0.5i) = -l, q(l + 0.5j)=j (30) 

The 4-soliton is generated by solving the Riemann-Hilbert 
linear system of equations with zero continuous spectrum [T] 
and can be seen in Fig. Eta). Figs. HJ a )-(d) show the discrete 
spectrum of the signal using various matrix-based methods. 
The relative accuracy of these schemes in estimating the 




N A 



(a) (b) 

Fig. 6. (a) Convergence of the discrete spectral amplitude for the rectangular 
pulse q(t) = 2rect(t) as a function of the number of sample points JV. Factor 
—j is not shown in the figure, (b) Continuous spectrum. 



eigenvalue A = 1 + 0.5j is shown in Fig. 13b). A very similar 
graph is obtained for other eigenvalues. Iterative methods fare 
similarly and their performance is shown in Fig. 0b). 

The convergence of the discrete spectral amplitudes q(Xj) is 
not quite satisfactory. Discrete spectral amplitudes associated 
with eigenvalues with small |5(Aj)| can be obtained with 
reasonable accuracy, although the convergence of q(Xj) is 
slower than the convergence of the eigenvalues themselves. 
On the other hand, discrete spectral amplitudes associated 
with eigenvalues with large |3(A)| are extremely sensitive 
to the location of eigenvalues and even slight changes in 
eigenvalues lead to radically different estimates for the spectral 
amplitudes. In fact, as the energy of the pulse is increased 
by having eigenvalues with large |5(A)|, the Riemann-Hilbert 
system becomes ill-conditioned. Therefore the discrete spectral 
amplitudes cannot generally be obtained using the methods 
discussed in this paper. It is illustrative to see the surface of 
|a(A)| in Fig. [21] The eigenvalues sometimes correspond to 
deep and narrow wells in the surface of |a(A)|, and sometimes 
they correspond to flat minima. In cases that they correspond 
to narrow wells, the derivative a' (A) is sensitive to the location 
of eigenvalues, leading to sensitivities in q(Xj). It is also clear 
from (fT~8b that a' (A) is proportional to A 2 and thus is sensitive 
to A. 
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Fig. 9. (a) Error in estimating the eigenvalue A = —1 + 0.25j in a 
4-soliton using search-based methods, (b) Error in estimating the discrete 
spectral amplitude \q\ = 1. 



Note that q(Xj) do not appear in the trace formula, and 
in particular they do not contribute to the signal energy. This 
part of the NFT controls the time center of the pulse and 
influences the signal phase too. Due to dependency on the time 
center of the signal and the fact that time center can hardly be 
used for digital transmission, the values of q(Xj) appear to be 
numerically chaotic and cannot carry much information. For 
this reason, we do not discuss these quantities in detail. 

VI. Nonlinear Fourier Transform of Pulses in 
Data Communications 

In this section, we use the numerical methods discussed 
in Section [TT] to compute the nonlinear Fourier transform 
of signals typically used in optical fiber transmission. The 
emphasis is on sine functions as they constitute signal degrees 
of freedom, but we also consider raised-cosine functions, sech 
signals, and Gaussian pulses. In particular, we study the effect 
of the amplitude and phase modulation on the structure of 
the nonlinear spectra. We will also discuss the spectrum of 
wavetrains formed by sine functions. 

Since the layer-peeling and the spectral methods give accu- 
rate results in estimating the nonlinear spectra, they are chosen 
for subsequent simulations. 

A. Amplitude and Phase Modulation of Sine Functions 

Fig. [TOl shows the spectrum of y(t) = Asinc(2t) under the 
amplitude modulation. It can be seen that a sine function is 
all dispersive as A is increased from zero, until about A = tt 
(IklUi = 1-275271") where a new eigenvalue emerges from 
the origin. Starting from A = 0, the continuous spectrum 
is a rectangle, resembling the ordinary Fourier transform 
— J-(q*(t))(2X). As A is increased, the continuous spectral 
function is narrowed until A = tt, where it looks like a 
delta function and its energy starts to deviate from the energy 
of the time domain signal. As A > tt is further increased, 
the dominant eigenvalue on the joj axis moves up until 
A = 1.277T, where Ai = 1.4234j and a new pair of eigenvalues 
emerges, starting from A23 = ±3.2 + 0.05j. When the newly 
created eigenvalues are not pronounced enough, for instance 
in this example when transiting from A = tt to A = 1.27 'tt, 
numerical algorithms have difficulties in determining whether 
these small emerging eigenvalues are part of the spectrum or 
not. Here it appears that for tt < A < 1.27 it there is just 
one dominant purely imaginary eigenvalue moving upward. 




Fig. 

A = 



10. Nonlinear Fourier transform of a sine function with amplitude 

1,2,3,4. 
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Fig. 11. Locus of eigenvalues of the sine function under amplitude 
modulation: (a) A = to A = 5, (b) A = to A = 20. 



At A = 1.27-7T, A23 emerge and move up in the complex 
plane as A is increased. An important observation is that 
the sine function appears to have not only purely imaginary 
eigenvalues, but also a pair of symmetric eigenvalues with 
nonzero real part emerging at high values of A; see Fig.fTTTb). 
This means that, for example, a sine function (viewed in the 
time domain) contains a stationary "central component" plus 
two small "side components" which travel to the left and right 
if the sine function is subject to the NLS flow. The locus of the 
eigenvalues of the function Asinc(2t) as a result of variations 
in A is given in Fig. [TT] 

It follows that the sine function is a simple example of a real 
symmetric pulse whose eigenvalues are not necessarily purely 
imaginary, as conjectured for a long time. However if q(t) is 
real, non-negative, and "single-lobe", then there are exactly 
N = [5 + ^[j 1 * 1 — ej eigenvalues, all purely imaginary |[T6l . 



Under phase modulation, in the form of adding a constant 
phase term to the signal, the eigenvalues and the magnitude 
of the continuous spectrum remained unchanged. Vertical shift 
in the phase of the continuous spectrum as a result of phase 
modulation can be seen in Fig. [T_2] 

We may also examine the effect of time-dependent phase 
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Fig. 12. Phase of the continuous spectram of a sine function when: (a) 
A = 4 (b) A = 4j. 



changes. The effect of linear chirp, of the form exp(jwi), is 
shown in picture Fig. Q~3] Linear chirp results in just a shift of 
the discrete and continuous spectrum to the left or the right, 
depending on the sign of the chirp. 

It is interesting to observe the effect of a quadratic chirp. 
The locus of eigenvalues that result due to changes in the 
quadratic phase qexp(juit 2 ) has been studied in lfl6ll for 
Gaussian pulses. In our sine function example, in the case that 
there is one discrete eigenvalue in the chirp-free case (such as 
when A — 4), increasing ui will move the eigenvalue on the jui 
axis upward, but then the eigenvalue move down again and is 
absorbed in the real axis. Fig. Q3fc) shows the moment before 
the eigenvalue is absorbed into the real axis. Note that the 
eigenvalues off the jui axis are considered to be spurious; their 
number increases as the number of sample points is increased. 

A more interesting behavior is observed when A = 12. 
Here, there are two eigenvalues on the jui axis: Ai s» 
10.0484j, A 2 = 5.5515j, together with A 3 , 4 = ±3.1315 + 
0.7462j (FiglHd)). As ui is increased, Ai and A2 move down 
and a fifth eigenvalue A5 emerges from the real axis and moves 
upwards on the jui axis. Eventually, at about u> — 41.41, A2 
and A5 "collide" and move out of the jui axis to the left and 
right. If ui is further increased, A2 and A5 are absorbed into 
the real axis; see Fig. [14] 

Collision of eigenvalues also occurs with time dilation. 
Signal q(t) = sinc(at) has 3 eigenvalues on the jui axis for 
a = 0.1, plus two small eigenvalues on two sides of the jui 
axis (Fig [ToT b)). As a is increased, the smaller eigenvalue on 
the jui axis comes down and a new eigenvalue is generated 
at the origin, moving upward. These two eigenvalues collide 
at 0A2j (a = 0.1330) and are diverted to the first and 
second quadrant, and eventually absorbed in the real axis at 
about m = ±0.32 (w = 0.1990), Fig. H6tc)-(d)). As a is 
decreased, more eigenvalues appear on the jui axis and fewer 
on the real axis (Fig. fToT e)-(f)). Note that the eigenvalues 
are not necessarily on the jui axis. For example, the signal 
y = sinc(0.1370i) clearly has eigenvalues Ai = 0.8684j, 
A 2 = 0.5797.7, A 3 , 4 = ±0.1055 + 0.1210*. 

Fig. [15] shows the nonlinear spectrum of a sine pulse under 
a quadratic chirp modulation, given by Ae^ 1 sinc(2i), is 
broadened as ui varies. 

The effect of time dilation on the continuous spectrum can 
be seen in Fig (T7] It can be observed that increasing band- 
width a, will increase the continuous range of real nonlinear 
frequencies, leading to bandwidth expansion. 
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Fig. 16. Locus of eigenvalues of sinc(at) as the bandwidth varies (a) 
from a = 0.1 to a = 0.6. (Eigenvalues with small SA are not shown here.) 
(b) Eigenvalues for a = 0.1. (c) Eigenvalues before collision and (d) after 
collision, (e) Eigenvalues for a = 0.06 before collision and (f) for a = 0.065 
after collision. 




Fig. 17. Bandwidth expansion in sinc(at) for (a) a = 0.06, a = 0.1, 
a = 0.3 (b) a = 1, a = 2, a = 3. 



B. Sine wavetrains 

The nonlinear spectrum of a wavetrain can take on a 
complicated form, just like its ordinary Fourier transform 
counterpart. Eigenvalues of a two-symbol train, for instance, 
depend on the amplitude and phase of the two signals, and 
their separation distance. 

We first analyze the case in which there are only two sine 
functions located at the fixed Nyquist distance from each other, 
i.e., y(t) = aisinc(2i + |) + a2sinc(2£ —5). For ai = a 2 = 2 
the spectrum consists of a single eigenvalue A = 0.3676j and a 
number of spurious eigenvalues as shown in Fig. [18} a). As the 
phase of a 2 is increased from 8 = to 6 = w, the eigenvalues 
moves off the jui axis to the left and a new eigenvalue emerges 
from the real axis in the first quadrant. Eigenvalues at 8 = u 
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Fig. 13. (a) Amplitude of the continuous spectrum with no carrier, (b) Amplitude of the continuous spectrum with carrier frequency ui = 5. The phase graph 
is also shifted similarly with no other change (AA = 2.5). (c) Locus of the eigenvalues of a sine function with amplitude A = 8 as the carrier frequency 
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■ 0, 10 and 30. (b) Phase of the 



are ±1.3908 + 0.3287i. The resulting locus of eigenvalues is 
shown in Fig. [T8l b). Figs [T8l c)-(d) depict similar graphs when 
both a\ and change. 

Next we study the locus of the discrete spectrum as a 
function of pulse separation for fixed amplitudes. If the ampli- 
tude of the sine functions is increased sufficiently, eigenvalues 
appear off the real axis and form a locus as the distance 
between pulses varies. Fig. [T9la) shows the locus of eigen- 
values of y(t) = 4sinc(2i + r) + 4sinc(2i — r) as r changes 
between zero to 5. At r = 0, eigenvalues are Ai = 6j, 
A 2 , 3 = ±2.3618 + 0.6476i and A 4j5 = ±3.2429 + 0.0815i As 
the distance between pulses is increased, Ai rapidly decreases, 



and at about r = 0.25, where Ai = 4.3j, the eigenvalues 
with non-zero real parts are absorbed into the real axis at 
3?A = —3. As t is further increased, Ai decreases further, 
until t = 0.4 where Ai = 2.4j and two new eigenvalues 
emerge at locations 3?A = ±3.12 going up and towards the 
ju axis. These eigenvalues return, before reaching the ju) axis, 
to be absorbed into the real axis, while new eigenvalues are 
generated again from the real axis. At r = 0.7 eigenvalues are 
A = 2j, ±1 +j. At some point, the newly created eigenvalues 
are not absorbed into the real axis, but they reach the jui axis 
and collide. A collision occurs, for instance, at r = 1.05. 
One of these eigenvalues goes down to be absorbed into the 
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smc(2t + |) 
< 7T, (c) ai = 2, 



origin, and the other one, interestingly, goes up to be united 
with the maximum eigenvalue on the juj axis (i.e., to create 
one eigenvalue with multiplicity two). Increasing the distance 
further does not change the location of this eigenvalue, which 
from now is fixed at A = lAj, but just changes the pattern 
of lower level eigenvalues. The collision does not occur when 
the amplitudes of the signals are smaller; see Fig. [19] (b) for 
the locus of the eigenvalues when the amplitude of the two 
sine functions is 2. 
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Fig. 19. (a) The locus of the discrete spectrum of y(t) = 4sinc(2t + 
r) + 4sinc(2t — r) as a function of < r < 5. (b) The locus of the discrete 
spectrum of y(t) = 2sinc(2t+r) + 2sinc(2t — r) as afunction of < r < 5. 



For wavetrains with a larger number of signals, the number 
of eigenvalues increases proportionally. We generate these 
wavetrains randomly and examine the region to which the 
spectrum is confined. Fig. [20] shows the locus of the discrete 
spectrum of all sine wavetrains with 16 signals. All 16 signal 
degrees of freedom in the bandlimited signal are modulated 
here. The effect of the bandwidth constraint in the nonlinear 
spectral domain can seen in this picture. 

C. Preservation of the spectrum of the NLS equation 

It is crucial to ensure that the spectrum found by the 
numerical methods, such as those discussed in the previous 
sections, is in fact correct. While it proved difficult to do so 




Fig. 20. Effect of the bandwidth constraint on the location of the eigenvalues 
of a sine wavetrain containing 16 pulses having random amplitudes. 

TABLE I 
Fiber Parameters 



D 

7 



17 ps/ (nm — km) 
1.27 W^km" 1 



Dispersion parameter 
nonlinearity parameter 



consistently and efficiently, there are various tests to increase 
one's confidence in the truth of the output of the numerical 
methods. Taking the inverse nonlinear Fourier transform in the 
continuous-time domain and comparing the resulting function 
in time with the original signal is generally quite cumbersome 
and not always feasible. One quick test is to examine a time 
frequency identity, such as the trace formula for n = 1,2,3,... 
as used in this paper. The first few conserved quantities in 
this identity can be written explicitly. One should allow higher 
tolerance values in the trace formula for large n, as the discrete 
terms in this identity involve A™ and thus are increasingly 
more sensitive to the eigenvalues. Another test is to subject the 
signal to the flow of an integrable equation, such as the NLS 
equation, and check that the discrete spectrum is preserved 
and the spectral amplitudes are scaled appropriately according 
to that equation. In this section, we let the signal propagate 
according to the NLS equation and compare the spectra at 
z = and z — C for various C. 

Fig.l2"T1shows examples of the spectra of a number of pulses 
at z = and z = C evolving according to the NLS equation 
(Q]l. The distances mentioned in the graphs in km correspond 
to a standard optical fiber with parameters in Table H] Note 
that in all these examples discrete the spectrum is completely 
preserved, and the continuous spectral amplitudes undergo a 
phase change properly. Compared to Gaussian and raised- 
cosine examples, whose nonlinear Fourier transform can be 
found easily, the discrete spectrum of sine functions is much 
more challenging to find. This is because the non-dominant 
eigenvalues off the jui axis have small imaginary parts for 
typical parameters and are not sufficiently distinguished. They 
also have large real parts, increasing the search region. Sine 
functions are thus not the best examples to illustrate the 
application of the NFT in optical fibers. We studied these ideal 
pulses primarily because of their fundamental utility in digital 
communications. 

VII. Conclusions 

In this paper, we have suggested and compared a variety 
of numerical methods for the computation of the nonlinear 
Fourier transform of a signal defined on the entire real line. 
A straightforward finite-difference discretization, such as the 
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forward discretization, does not often produce satisfactory 
results. Among the methods studied in this paper, the layer- 
peeling and spectral methods gave accurate results in estimat- 
ing the continuous and discrete spectrum over a wide class of 
examples. 

Given a waveform without having prior knowledge of the 
location of the discrete eigenvalues, we suggest the use of 
matrix-based methods to compute the discrete spectrum. If, 
on the other hand, the location of the eigenvalues is known 
approximately (as in data-communication problems, where the 
eigenvalues are chosen at the transmitter from a finite set) a 
search-based method is recommended. 

Although the eigenvalues and the continuous spectral func- 
tion can be calculated with great accuracy, the discrete spectral 
amplitudes are quite sensitive to the location of the eigenval- 
ues, even in the absence of noise. These discrete amplitudes 
control the time center of the pulse, and are therefore sensitive 
to timing jitter. For data communication purposes it follows 
that, whereas the presence or absence of the eigenvalue itself 
may allow for robust information transmission, encoding in- 
formation in the time center of the pulse, i.e., in the discrete 
spectral amplitudes, is unlikely to be viable. 

Using these numerical methods, we studied the influence of 
various signal parameters on the nonlinear Fourier transform 
of a number of pulses commonly used in data communications. 
We found, for example, that the spectrum of an isolated nor- 
malized sine function with amplitude A is purely continuous 
for the A < ir. However, as the pulse amplitude is increased, 
dominant eigenvalues appear on the ju> axis, together with 
pairs of symmetric eigenvalues having nonzero real part. 

In general, amplitude variations result in variations in the 
location of the eigenvalues and the shape of the continuous 
spectrum. Eigenvalues follow particular trajectories in the 
complex plane. Phase variations, on the other hand, influ- 
ence only the phase of the spectrum, not the location of 
the eigenvalues. One important observation, which may be 
beneficial for the design of data communication systems, is 
that the nonlinear spectrum of bandlimited pulses appears to 
be confined to a vertical strip in the complex plane with a 
width proportion to the signal bandwidth. 

This paper has only scratched the surface of a poten- 
tially rich research area. The development of efficient and 
robust numerical techniques suitable for various engineering 
applications of the nonlinear Fourier transform will require 
significant additional effort. A problem of particular interest 
is the development of a "fast" nonlinear Fourier transform 
method that would be the analog of the FFT 
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